Notes

At first I used HTODemux for Day 3 and MULTIseqDemux for Day 15 because the ridge plots had better separation at Day 3 that way. However, we get 1664 more singlets with MULTIseqDemux at Day 3 and the hashtag assignment is almost entirely in agreement between the two tools. So I’ve decided to use MULTIseqDemux for both timepoint for consistency.

Load libraries and data

library(Seurat)
library(tidyverse)
library(knitr)
library(cowplot)
library(ggplot2)

set.seed(4)

script_output_dir <- file.path(here::here(), "output")

# Create folders for data
if(!dir.exists(file.path(script_output_dir))) {
  cat(sprintf("Creating folder %s\n", file.path(script_output_dir, "processed_data")))
  dir.create(file.path(script_output_dir, "processed_data"), recursive = T)
  cat(sprintf("Creating folder %s\n", file.path(script_output_dir, "plots")))
  dir.create(file.path(script_output_dir, "plots"), recursive = T)
}

Load data into Seurat objects

# Linux
day3_f <- "/media/emmabishop/5TBSharedStorage/project_data/2022_BCGChallenge/Round1/20230323_Day3_aggr/outs/count/filtered_feature_bc_matrix.h5"
day15_f <- "/media/emmabishop/5TBSharedStorage/project_data/2022_BCGChallenge/Round1/20230207_Day15/per_sample_outs/20230203_Seshadri_UW/count/sample_filtered_feature_bc_matrix.h5"

Make objects

make_seurat_obj <- function(infile, projectname) {
  dat <- Read10X_h5(infile)
  
  # Remove ".1" from gene names
  rownames(dat$`Antibody Capture`) <- rownames(dat$`Antibody Capture`) %>%
    str_replace(., fixed(".1"), "")
  
  # Extract matrices of ADT and HTO info 
  ADT <- !grepl("Hashtag", dat$`Antibody Capture`@Dimnames[[1]])
  ADT_counts <- dat$`Antibody Capture`[ADT,]
  HTO <- grepl("Hashtag", dat$`Antibody Capture`@Dimnames[[1]])
  HTO_counts <- dat$`Antibody Capture`[HTO,]
  
  # Make object
  obj <- CreateSeuratObject(counts = dat$`Gene Expression`, project = projectname, min.cells = 3)
  obj[["ADT"]] <-  CreateAssayObject(counts = ADT_counts)
  obj[["HTO"]] <-  CreateAssayObject(counts = HTO_counts)
  Assays(obj)
  
  return(obj)
}

day3 <- make_seurat_obj(day3_f, "Day3")
day15 <- make_seurat_obj(day15_f, "Day15")

# Starting counts
d3_start_count <- length(day3$orig.ident)
d15_start_count <- length(day15$orig.ident)

Raw HTO counts for each day

rowSums(day3$HTO) %>%
  as.data.frame() %>%
  kable(align = "l", caption = "Day 3")
Day 3
.
Hashtag1 805481
Hashtag2 1928906
Hashtag3 1490418
Hashtag4 1953263
Hashtag5 1019046
Hashtag6 3034373
Hashtag7 5928
Hashtag8 791903
Hashtag9 2986667
Hashtag10 2705637 
rowSums(day15$HTO) %>%
  as.data.frame() %>%
  kable(align = "l", caption = "Day 15")
Day 15
.
Hashtag1 4067421
Hashtag3 1405752
Hashtag4 204878
Hashtag5 3176998
Hashtag6 1102579
Hashtag7 8374125
Hashtag9 3458347
Hashtag10 2079834

Normalize data

day3_norm <- NormalizeData(day3, assay = "HTO", normalization.method = "CLR", margin = 2)
## Normalizing across cells
day15_norm <- NormalizeData(day15, assay = "HTO", normalization.method = "CLR", margin = 2)
## Normalizing across cells

HTO distributions for each day

hto_distr <- function(inobj) {
  HTO_counts <- inobj$HTO@data
  to_plot <- HTO_counts %>%
    as.data.frame() %>%
    rownames_to_column("hashtag") %>%
    pivot_longer(!hashtag) %>%
    select(-c(name))
  
  h <- ggplot(to_plot, aes(value)) +
    geom_histogram() +
    theme_classic() +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    facet_wrap(~hashtag, scales = "free")
  
  return(h)
}

Day 3

d3_hto <- hto_distr(day3_norm)
d3_hto + ggtitle("Day 3")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Day 15

d15_hto <- hto_distr(day15_norm)
d15_hto + ggtitle("Day 15")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Demultiplex

demux_multiseq <- function(norm_obj) {
  # Demultiplex with MULTIseqDemux 
  dm <- MULTIseqDemux(norm_obj, assay = "HTO", autoThresh = TRUE)
  i1 <- grepl("Hashtag", dm@meta.data$MULTI_ID)
  i2 <- grepl("Doublet", dm@meta.data$MULTI_ID)
  i3 <- grepl("Negative", dm@meta.data$MULTI_ID)
  # Get overall Singlet, Doublet, Negative classification
  dm@meta.data$MULTI_ID_global <- NULL
  dm@meta.data$MULTI_ID_global[i1] <- "Singlet"
  dm@meta.data$MULTI_ID_global[i2] <- "Doublet"
  dm@meta.data$MULTI_ID_global[i3] <- "Negative"
  dm@meta.data$MULTI_ID_global <- factor(dm@meta.data$MULTI_ID_global, 
                                         levels = c("Doublet", 
                                                    "Negative", 
                                                    "Singlet"))
  return(dm)
}


day3_dm <- demux_multiseq(day3_norm)
day3_dm@meta.data$MULTI_ID <- factor(day3_dm@meta.data$MULTI_ID,
                                levels = c("Hashtag1", "Hashtag2", "Hashtag3", 
                                           "Hashtag4", "Hashtag5", "Hashtag6", 
                                           "Hashtag7", "Hashtag8", "Hashtag9", 
                                           "Hashtag10", "Doublet", "Negative"))

day15_dm <- demux_multiseq(day15_norm)
day15_dm@meta.data$MULTI_ID <- factor(day15_dm@meta.data$MULTI_ID,
                                levels = c("Hashtag1", "Hashtag3", 
                                           "Hashtag4", "Hashtag5", "Hashtag6", 
                                           "Hashtag7", "Hashtag9", 
                                           "Hashtag10", "Doublet", "Negative"))

Day 3

Missing hashtags:

  • Hashtag 7 (ptid 12) - why?
table(day3_dm$MULTI_ID_global)
## 
##  Doublet Negative  Singlet 
##     1583     1828    13053
table(day3_dm$MULTI_ID)
## 
##  Hashtag1  Hashtag2  Hashtag3  Hashtag4  Hashtag5  Hashtag6  Hashtag7  Hashtag8 
##       762      1250      1257      1612      1423      1404         5      1461 
##  Hashtag9 Hashtag10   Doublet  Negative 
##      1634      2245      1583      1828

Ridge plot using MULTI_ID

Idents(day3_dm) <- "MULTI_ID"
RidgePlot(day3_dm, assay = "HTO", ncol = 3,
          features = c("Hashtag1", "Hashtag2", "Hashtag3", "Hashtag4", "Hashtag5",
                       "Hashtag6", "Hashtag8", "Hashtag9", "Hashtag10"))

Day 15

Missing hashtags:

  • Hashtag 2 (ptid 5) - Expected
  • Hashtag 3 (ptid 7) - why?
  • Hashtag 4 (ptid 8) - why?
  • Hashtag 8 (ptid 12) - Expected
table(day15_dm$MULTI_ID_global)
## 
##  Doublet Negative  Singlet 
##      466      484     5449
table(day15_dm$MULTI_ID)
## 
##  Hashtag1  Hashtag3  Hashtag4  Hashtag5  Hashtag6  Hashtag7  Hashtag9 Hashtag10 
##      1191        16         9      1286       300      1442       730       475 
##   Doublet  Negative 
##       466       484

Ridge plot using MULTI_ID

Idents(day15_dm) <- "MULTI_ID"
RidgePlot(day15_dm, assay = "HTO", ncol = 3,
          features = c("Hashtag1", "Hashtag3", "Hashtag4", "Hashtag5",
                       "Hashtag6", "Hashtag7", "Hashtag9", "Hashtag10"))

Label by HTO

label <- function(dm) {
  i1 <- grepl("^Hashtag1$", dm@meta.data$MULTI_ID)
  i2 <- grepl("Hashtag2$", dm@meta.data$MULTI_ID)
  i3 <- grepl("^Hashtag3$", dm@meta.data$MULTI_ID)
  i4 <- grepl("^Hashtag4$", dm@meta.data$MULTI_ID)
  i5 <- grepl("^Hashtag5$", dm@meta.data$MULTI_ID)
  i6 <- grepl("^Hashtag6$", dm@meta.data$MULTI_ID)
  i7 <- grepl("^Hashtag7$", dm@meta.data$MULTI_ID)
  i8 <- grepl("^Hashtag8$", dm@meta.data$MULTI_ID)
  i9 <- grepl("^Hashtag9$", dm@meta.data$MULTI_ID)
  i10 <- grepl("^Hashtag10$", dm@meta.data$MULTI_ID)
  i11 <- grepl("Doublet", dm@meta.data$MULTI_ID)
  i12 <- grepl("Negative", dm@meta.data$MULTI_ID)

  # PTID number
  dm@meta.data$PTID <- NULL
  dm@meta.data$PTID[i1] <- "1"
  dm@meta.data$PTID[i2] <- "5"
  dm@meta.data$PTID[i3] <- "7"
  dm@meta.data$PTID[i4] <- "8"
  dm@meta.data$PTID[i5] <- "9"
  dm@meta.data$PTID[i6] <- "10"
  dm@meta.data$PTID[i7] <- "11"
  dm@meta.data$PTID[i8] <- "12"
  dm@meta.data$PTID[i9] <- "13"
  dm@meta.data$PTID[i10] <- "16"
  dm@meta.data$PTID[i11] <- "Doublet"
  dm@meta.data$PTID[i12] <- "Negative"
  dm@meta.data$PTID <- factor(dm@meta.data$PTID, 
                              levels = c("1", "5", "7", "8", "9", "10", "11",
                                         "12", "13", "16", "Doublet", "Negative"))
  # PTID sex
  dm@meta.data$SEX <- NULL
  dm@meta.data$SEX[i1] <- "M"
  dm@meta.data$SEX[i2] <- "F"
  dm@meta.data$SEX[i3] <- "I"
  dm@meta.data$SEX[i4] <- "M"
  dm@meta.data$SEX[i5] <- "F"
  dm@meta.data$SEX[i6] <- "M"
  dm@meta.data$SEX[i7] <- "M"
  dm@meta.data$SEX[i8] <- "M"
  dm@meta.data$SEX[i9] <- "M"
  dm@meta.data$SEX[i10] <- "F"
  dm@meta.data$SEX[i11] <- "Doublet"
  dm@meta.data$SEX[i12] <- "Negative"

  # PTID age
  dm@meta.data$AGE <- NULL
  dm@meta.data$AGE[i1] <- 29
  dm@meta.data$AGE[i2] <- 40
  dm@meta.data$AGE[i3] <- 40
  dm@meta.data$AGE[i4] <- 41
  dm@meta.data$AGE[i5] <- 36
  dm@meta.data$AGE[i6] <- 28
  dm@meta.data$AGE[i7] <- 36
  dm@meta.data$AGE[i8] <- 43
  dm@meta.data$AGE[i9] <- 28
  dm@meta.data$AGE[i10] <- 20
  dm@meta.data$AGE[i11] <- 0
  dm@meta.data$AGE[i12] <- 0

  # Arm metadata
  dm@meta.data$ARM <- NULL
  dm@meta.data$ARM[i1] <- "Non-INH"
  dm@meta.data$ARM[i2] <- "Non-INH"
  dm@meta.data$ARM[i3] <- "Non-INH"
  dm@meta.data$ARM[i4] <- "Non-INH"
  dm@meta.data$ARM[i5] <- "INH"
  dm@meta.data$ARM[i6] <- "INH"
  dm@meta.data$ARM[i7] <- "INH"
  dm@meta.data$ARM[i8] <- "Non-INH"
  dm@meta.data$ARM[i9] <- "INH"
  dm@meta.data$ARM[i10] <- "INH"
  dm@meta.data$ARM[i11] <- "Doublet"
  dm@meta.data$ARM[i12] <- "Negative"
  dm@meta.data$ARM <- factor(dm@meta.data$ARM, levels = c("Non-INH", "INH"))
  
  return(dm)
}

day3_labelled <- label(day3_dm)
day15_labelled <- label(day15_dm)

saveRDS(day3_labelled, file = file.path(script_output_dir, "processed_data/1_dm_all_d3.rds"))
saveRDS(day15_labelled, file = file.path(script_output_dir, "processed_data/1_dm_all_d15.rds"))

QC

ncount RNA (singlets, doublets, negative)

Idents(day3_labelled) <- "MULTI_ID_global"
v1 <- VlnPlot(day3_labelled, features = "nCount_RNA", pt.size = 0.1, log = TRUE) +
  ylim(0, 40000) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 3 nCount_RNA") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
Idents(day15_labelled) <- "MULTI_ID_global"
v2 <- VlnPlot(day15_labelled, features = "nCount_RNA", pt.size = 0.1, log = TRUE) +
  ylim(0, 40000) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 15 nCount_RNA") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
cowplot::plot_grid(v1, v2)
## Warning: Removed 198 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 198 rows containing missing values (`geom_point()`).

ncount HTO (singlets, doublets, negative)

v3 <- VlnPlot(day3_labelled, features = "nCount_HTO", pt.size = 0.1, log = TRUE) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 3 nCount_HTO") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
v4 <- VlnPlot(day15_labelled, features = "nCount_HTO", pt.size = 0.1, log = TRUE) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 15 nCount_HTO") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
cowplot::plot_grid(v3, v4)

ncount ADT (singlets, doublets, negative)

v5 <- VlnPlot(day3_labelled, features = "nCount_ADT", pt.size = 0.1, log = TRUE) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 3 nCount_ADT") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
v6 <- VlnPlot(day15_labelled, features = "nCount_ADT", pt.size = 0.1, log = TRUE) +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Day 15 nCount_ADT") +
  NoLegend()
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
cowplot::plot_grid(v5, v6)

Total Cells per PTID and arm

Get singlets

Idents(day3_labelled) <- "MULTI_ID_global"
Idents(day15_labelled) <- "MULTI_ID_global"
d3_singlet <- subset(day3_labelled, idents = "Singlet")
d15_singlet <- subset(day15_labelled, idents = "Singlet")

# Get singlet by HTO counts
d3_hto_singlet_count <- length(d3_singlet$orig.ident)
d15_hto_singlet_count <- length(d15_singlet$orig.ident)

# Start QC barplot dataframe
d3_qc_counts <- data.frame(d3_start_count, d3_hto_singlet_count)
d15_qc_counts <- data.frame(d15_start_count, d15_hto_singlet_count)

# These are necessary for the next script
saveRDS(d3_qc_counts, file = file.path(script_output_dir, "processed_data/2_d3_qc_df.rds"))
saveRDS(d15_qc_counts, file = file.path(script_output_dir, "processed_data/2_d15_qc_df.rds"))
saveRDS(d3_singlet, file = file.path(script_output_dir, "processed_data/2_dm_singlets_d3.rds"))
saveRDS(d15_singlet, file = file.path(script_output_dir, "processed_data/2_dm_singlets_d15.rds"))
df_day3 <- table(d3_singlet$orig.ident, d3_singlet$PTID) %>%
  as.data.frame()

df_out <- table(d15_singlet$orig.ident, d15_singlet$PTID) %>%
  as.data.frame() %>%
  bind_rows(df_day3) %>%
  filter(Var2 != "Doublet" & Var2 != "Negative") %>%
  dplyr::rename("PTID" = "Var2") %>%
  dplyr::rename("Timepoint" = "Var1") %>%
  dplyr::rename("# Cells" = "Freq") %>%
  mutate(Timepoint = as.factor(Timepoint))
df_out$Timepoint <- factor(df_out$Timepoint, levels = c("Day3", "Day15"))

# Barplot
ggplot(df_out, aes(x = PTID, y = `# Cells`, fill = Timepoint)) +
  geom_bar(position = "dodge", stat = "identity", width = 0.7, color="black") +
  theme_bw() +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  ggtitle("Number of Singlet Cells per PTID")

Number of Non-INH vs INH ptids and cells at each day

dfout2 <- df_out %>%
  mutate(Arm = case_when(
    PTID %in% c("1", "5", "7", "8", "12") ~ "Non-INH",
    .default = "INH"
  )) %>%
  group_by(Timepoint, Arm) %>%
  mutate(`Total Cells` = sum(`# Cells`)) %>%
  mutate(`# PTIDs with > 20 cells` = case_when(
    Timepoint == "Day15" & Arm == "Non-INH" ~ "n = 1",
    Timepoint == "Day15" & Arm == "INH" ~ "n = 5",
    Timepoint == "Day3" & Arm == "Non-INH" ~ "n = 5",
    .default = "n = 4"
  )) %>%
  ungroup()

ggplot(dfout2, aes(x = Timepoint, y = `Total Cells`, fill = Arm)) +
  geom_bar(position = "dodge", stat = "identity", width = 0.7) +
  theme_bw() +
  theme(text = element_text(size = 16),
        axis.text = element_text(size = 16)) +
  scale_fill_manual(values = c("orange", "blue4")) +
  geom_text(aes(label=`# PTIDs with > 20 cells`), 
            position=position_dodge(width = 0.7), 
            vjust=-0.25) +
  ggtitle("Total Singlets per Arm (n is # PTIDs with > 20 cells)")

dfout2 %>%
  group_by(Timepoint, Arm, `# Cells`, PTID) %>%
  tally() %>%
  select(-c(n)) %>%
  kable(align = "l", caption = "Singlets per Timepoint, Arm, PTID")
Singlets per Timepoint, Arm, PTID
Timepoint Arm # Cells PTID
Day3 INH 5 11
Day3 INH 1404 10
Day3 INH 1423 9
Day3 INH 1634 13
Day3 INH 2245 16
Day3 Non-INH 762 1
Day3 Non-INH 1250 5
Day3 Non-INH 1257 7
Day3 Non-INH 1461 12
Day3 Non-INH 1612 8
Day15 INH 300 10
Day15 INH 475 16
Day15 INH 730 13
Day15 INH 1286 9
Day15 INH 1442 11
Day15 Non-INH 0 5
Day15 Non-INH 0 12
Day15 Non-INH 9 8
Day15 Non-INH 16 7
Day15 Non-INH 1191 1 
dfout2 %>%
  group_by(Timepoint, Arm, `Total Cells`, `# PTIDs with > 20 cells`) %>%
  tally() %>%
  select(-c(n)) %>%
  kable(align = "l", caption = "Singlets per Timepoint, Arm")
Singlets per Timepoint, Arm
Timepoint Arm Total Cells # PTIDs with > 20 cells
Day3 INH 6711 n = 4
Day3 Non-INH 6342 n = 5
Day15 INH 4233 n = 5
Day15 Non-INH 1216 n = 1

tSNE plots

Get Singlets and Doublets

day3_subset <- subset(day3_labelled, idents = "Negative", invert = TRUE)
day15_subset <- subset(day15_labelled, idents = "Negative", invert = TRUE)

tSNE embedding of the HTO data

tsne_hto <- function(dm_subset) {
  DefaultAssay(dm_subset) <- "HTO"
  dm_subset <- ScaleData(dm_subset, features = rownames(dm_subset),
                                   verbose = FALSE)
  dm_subset <- RunPCA(dm_subset, features = rownames(dm_subset), approx = FALSE)
  dm_subset <- RunTSNE(dm_subset, dims = 1:8, perplexity = 100)
  DimPlot(dm_subset)
  
  return(dm_subset)
}

day3_tsne <- tsne_hto(day3_subset)
day15_tsne <- tsne_hto(day15_subset)

Doublets vs Singlets

p1 <- DimPlot(day3_tsne) + ggtitle("Day 3")
p2 <- DimPlot(day15_tsne) + ggtitle("Day 15")

p1 + p2

Singlets by PTID

ptid_colors <- c("1" = "salmon", "5" = "orange", "7" = "#7570B3",
                 "8" = "#0CB702", "9" = "#A6761D", "10" = "#00BFC4",
                 "11" = "#C77CFF", "12" = "#FF61CC", "13" = "darkgreen",
                 "16" = "darkblue")

# Get tSNE singlets
d3_tsne_singlet <- subset(day3_tsne, idents = "Singlet")
d15_tsne_singlet <- subset(day15_tsne, idents = "Singlet")

Idents(d3_tsne_singlet) <- "PTID"
Idents(d15_tsne_singlet) <- "PTID"

p3 <- DimPlot(d3_tsne_singlet, cols = ptid_colors) + ggtitle("Day 3")
p4 <- DimPlot(d15_tsne_singlet, cols = ptid_colors) + ggtitle("Day 15")

p3 + p4

Session Info

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Los_Angeles
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.1      knitr_1.45         lubridate_1.9.3    forcats_1.0.0     
##  [5] stringr_1.5.0      dplyr_1.1.3        purrr_1.0.2        readr_2.1.5       
##  [9] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.4.4      tidyverse_2.0.0   
## [13] Seurat_5.0.1       SeuratObject_5.0.1 sp_2.1-1          
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.15.0      jsonlite_1.8.7        
##   [4] magrittr_2.0.3         ggbeeswarm_0.7.2       spatstat.utils_3.0-4  
##   [7] farver_2.1.1           rmarkdown_2.25         vctrs_0.6.4           
##  [10] ROCR_1.0-11            spatstat.explore_3.2-5 htmltools_0.5.7       
##  [13] sass_0.4.7             sctransform_0.4.1      parallelly_1.36.0     
##  [16] KernSmooth_2.23-22     bslib_0.5.1            htmlwidgets_1.6.2     
##  [19] ica_1.0-3              plyr_1.8.9             plotly_4.10.3         
##  [22] zoo_1.8-12             cachem_1.0.8           igraph_1.5.1          
##  [25] mime_0.12              lifecycle_1.0.4        pkgconfig_2.0.3       
##  [28] Matrix_1.6-4           R6_2.5.1               fastmap_1.1.1         
##  [31] fitdistrplus_1.1-11    future_1.33.0          shiny_1.7.5.1         
##  [34] digest_0.6.33          colorspace_2.1-0       patchwork_1.1.3       
##  [37] rprojroot_2.0.3        tensor_1.5             RSpectra_0.16-1       
##  [40] irlba_2.3.5.1          labeling_0.4.3         progressr_0.14.0      
##  [43] fansi_1.0.6            spatstat.sparse_3.0-3  timechange_0.3.0      
##  [46] httr_1.4.7             polyclip_1.10-6        abind_1.4-5           
##  [49] compiler_4.3.3         here_1.0.1             bit64_4.0.5           
##  [52] withr_3.0.0            fastDummies_1.7.3      highr_0.10            
##  [55] MASS_7.3-60.0.1        tools_4.3.3            vipor_0.4.5           
##  [58] lmtest_0.9-40          beeswarm_0.4.0         httpuv_1.6.12         
##  [61] future.apply_1.11.0    goftest_1.2-3          glue_1.7.0            
##  [64] nlme_3.1-163           promises_1.2.1         grid_4.3.3            
##  [67] Rtsne_0.16             cluster_2.1.6          reshape2_1.4.4        
##  [70] generics_0.1.3         hdf5r_1.3.8            gtable_0.3.4          
##  [73] spatstat.data_3.0-3    tzdb_0.4.0             data.table_1.15.0     
##  [76] hms_1.1.3              utf8_1.2.4             spatstat.geom_3.2-7   
##  [79] RcppAnnoy_0.0.21       ggrepel_0.9.4          RANN_2.6.1            
##  [82] pillar_1.9.0           spam_2.10-0            RcppHNSW_0.5.0        
##  [85] later_1.3.1            splines_4.3.3          lattice_0.22-5        
##  [88] bit_4.0.5              survival_3.5-8         deldir_1.0-9          
##  [91] tidyselect_1.2.0       miniUI_0.1.1.1         pbapply_1.7-2         
##  [94] gridExtra_2.3          scattermore_1.2        xfun_0.40             
##  [97] matrixStats_1.0.0      stringi_1.8.3          lazyeval_0.2.2        
## [100] yaml_2.3.8             evaluate_0.22          codetools_0.2-19      
## [103] cli_3.6.2              uwot_0.1.16            xtable_1.8-4          
## [106] reticulate_1.34.0      munsell_0.5.0          jquerylib_0.1.4       
## [109] Rcpp_1.0.11            globals_0.16.2         spatstat.random_3.2-1 
## [112] png_0.1-8              ggrastr_1.0.2          parallel_4.3.3        
## [115] ellipsis_0.3.2         dotCall64_1.1-0        listenv_0.9.0         
## [118] viridisLite_0.4.2      scales_1.3.0           ggridges_0.5.4        
## [121] leiden_0.4.3           rlang_1.1.1